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Abstract 

In this paper we employ two distinct approaches - all-electron ab initio method and the spherical 

^jellium background model- within time dependent density functional theory to calculate the long 
-I— • 

^ ^ange dipole-dipole dispersion coefficient (van der Waals coefficient) Cg of sodium atom clusters 
• rH containing even number of atoms ranging from 2 to 20 atoms. The dispersion coefficients are 
I^Mained .a Ca.ni.PoM. .Mion. A,> ..e »,c.a«o. a. ca„ie. ... ,oca, de„* ap- 
L^j|)roximation for exchange-correlation energy functional. These two sets of results are compared to 
assess the accuracy of jellium based results and to ascertain the effect of detail ionic structure of 
3^ the clusters on the van der Waals coefficient. 
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I. INTRODUCTION 



Contribution of long-range dispersive van der Waals forces to the interaction between 
two molecules is quite significant. These forces play an important role in the description 
of many physical and chemical phenomena. Some of these are adhesion, surface tension, 
physical adsorption chemistry of rare gases Q|, molecular chemistry in the interstellar 



medium [s] etc.. The correlation between the electron density fluctuation at widely separated 
locations give rise to these forces. The van Waals interaction between two molecules A and 
B in the lowest order is represented by the orient at ionally averaged potential 

VABiR) = -§ (1) 

where R is the intermolecular distance and Cq is the orientationally averaged van der Waals 
coefficient. The coefficient Ce describes the dipole-dipole interactions between the two molec- 
ular species. 

Metal clusters, specifically those of alkali atoms Na and K, played an important role in the 
development of cluster physics as a branch of modern physics and chemistry. The knowledge 
of their van der Waals coefficient Cq is useful for the description of cluster-cluster collisions 
4] and also for characterizing the orientation of clusters in the bulk matter [sl, Q]. However, 
only very few papers devoted to the calculation of van der Waals coefficients exist in the 
literature. In Refs. 0,0], time dependent Kohn-Sham (TDKS) equation of time dependent 
density functional theory (TDDFT) within the spherical jellium background model (SJBM) 
was employed to calculate the van der Waals coefficients. On the other hand, in Ref. a 
purely particle- and current-density based hydrodynamic formalism of TDDFT within SJBM 
has been applied to calculate the coefficients. This model is well suited for the description of 
the alkali metal clusters (Li, Na, and K) and correctly explains the greater stability of these 
clusters with magic number 2, 8, 20, 34, 40, •■■ of atoms due to closing of the electronic shell 
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Although the hydrodynamic approach is less accurate than the TDKS approach 
due to the approximate nature of the kinetic energy functional employed in it, nonetheless 
it is computationally advantageous for larger clusters in comparison to the orbital-based 
TDKS approach. 

The SJBM replaces the discrete ionic structure of clusters by a spherically symmetric 
uniform positive charge background and thus making it possible to carry out calculations 
for the optical response properties of reasonably large clusters of around 100 atoms within the 



TDKS approach 12|. In contrast to this the hydrodynamic approach in conjunction with 
SJBM allows calculation of optical response properties of very large alkali atom clusters 



containing 10,000 atoms without much difficulty [9|, [l3|, [ij, [l5[. In past ten years or so 



several all-electron ab initio calculations devoted to the ground state and the optical response 



properties of sodium clusters taking into account the actua 
sodium atoms have also been reported in the literature 
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26| . However, these calculations are restricted to small sized clusters containing maximum 



up to 20 sodium atoms. Moreover, these calculations are computationally more expensive 
and require substantially more resources than the jellium based calculations. It is only very 
recently that the all-electron ab initio calculations of the van der Waals coefficient Cg of 
small sized closed shell sodium clusters containing up to 20 atoms have been reported in 
the literature 



271, 



281]. The main aim of this paper is to make a systematic comparison of 



the values of Cg obtained by employing the jellium based model and all-electron ab initio 
method within the TDDFT. This comparison will enable us to assess the accuracy of jellium 
model in predicting the values of van der Waals coefficient and such a study is important, 
as the jellium model often turns out to be much more efficient to handle particularly when 
dealing with larger cluster systems. 

Before proceeding further, it is important to note that density functional theory (DFT) 
in principle should give the exact ground-state properties including the long range van der 
Waals energies. However, the widely used local density approximation (LDA) and general- 



ized gradient approximations (GGA) [29 



30 



31| exchange-correlation (XC) functionals fail 



to reproduce the van der Waals energies. This is because the LDA and GGA XC functionals 
cannot completely simulate the correlated motion of electrons arising from Coulomb inter- 
action between distant non overlapping electronic systems. It is only recently that attempts 
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34| have been made to obtain van der Waals energies directly from the ground-state 



energy functional by correcting the long range nature of the effective Kohn-Sham potential. 
On the other hand, it is possible to make reliable estimates of the van der Waals coefficient 
Cq by using expressions which relate this coefficient to the frequency dependent dipole po- 
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361]. We follow the latter route for the calculation 



larizabilities at imaginary frequencies 
of these coefficients. 

The paper is organized as follows: In section II we discuss the theoretical method and 
the expressions employed to calculate the van der Waals coefficient Cq from the frequency 



3 



dependent dipole polarizability. Results of our calculations are presented in Section III. The 
paper is concluded in Section IV. 



II. METHOD OF CALCULATION 



In order to calculate the van der Waals coefficient Ce, we make use of the Casimir- Polder 
expression which relates Cq to the frequency dependent dipole polarizability evaluated at 
imaginary frequency. In accordance with this expression the orientation averaged dispersion 
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coefficient between two molecules A and B is given by 

3 f°° 
Jo 



36| 



(2) 



where <yj{iuj) is the isotropic average dipole polarizability of the j-th molecule and is given 
by 

aj{uj) = ^- . (3) 



In the above expression a^xi^), ctyyi^^) and 0:2^(0;) are diagonal elements of the dipole polar- 
izability tensor. Therefore, the calculation of dispersion coefficient Cg involves determining 
frequency dependent dipole polarizability tensor followed by the evaluation of the quadra- 
ture. For the determination of the frequency dependent polarizability, we use two different 
methods (i) the hydrodynamic formalism of TDDFT in conjunction with the SJBM and (ii) 
all-electron ab initio method also in the realm of TDDFT. For details, we refer the reader 



to Refs. , [15(1 for the hydrodynamic formalism and Ref. [37|] for the ab initio method 

of calculation. In the following, we describe the choice of parameters made in this paper to 
perform the calculations of the dispersion coefficient Cq by two methods. 



n 



A. Hydrodynamic-SJBM Formalism 

The basic idea of the SJBM is to replace the distribution of ionic cores by a constant 
spherically symmetric positive background of jellium density ?T-+(r) in a finite volume. This 
positive background provides the attractive potential for the valence electrons contributed 
by each atom of the cluster. The ground-state electronic distribution is then obtained 
by employing either orbital based Kohn-Sham or purely density based extended Thomas- 
Fermi (ETF) formalisms of DFT. The calculation of the frequency dependent polarizabilty 
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within the SJBM can be accomphshed by TDDFT. Like its time independent counterpart 
TDDFT can be formulated either in terms of the time dependent orbitals (TDKS formahsm) 
or in terms of the particle and the current densities (hydrodynamic formalism). In this 
paper, we employ the hydrodynamic approach formalism of TDDFT in conjunction with 
the SJBM to calculate the frequency dependent polarizabilities at imaginary frequencies. 
The basic dynamical variables of the hydrodynamic theory are the time-dependent density 
p(r, t) and the velocity potential S{r,t). The velocity of the electron fluid is given by 
v(r, t) = —VS{r,t). Thus the total time-averaged energy can be expressed in terms of 
these two variables. For our purpose we need to evaluate the second-order change in the 
time-averaged energy as this is related directly to the frequency dependent average dipole 
polarizability by the relation 

a{uj) = -AE^^\iu) (4) 

where u is the frequency of the applied electromagnetic field. It has been shown that 
the second-order change in the time-average energy E'^'^\uj) is stationary with respect to 
the variations in the first-order induced density p^^\y^uj) and the induced current-density 
S^^\r,uj) Q. Consequently, E^^) {uo) can be determined by choosing appropriate variational 
forms for /3^^-'(r, uj) and S'^^\y, uj) and making i?*^^) stationary with respect to the parameters 
of p''^'*(r, uj) and 5'''^'*(r, oj). For our calculations we use ten parameters each for p'^'^^v^ uj) and 
S'^^\y,uj) and check the convergence by adding more parameters. The variational method 
is also applicable for the imaginary frequencies 9|] as it requires replacing tu^ by —uj"^ in the 
expression of E^'^\uj). This allows us to determine dynamic multipolarizability at imaginary 
frequencies {a{iuj)) by exactly the same procedure as employed for getting a{uj). 

It is well known that a TDDFT based response property calculation requires approximat- 
ing the XC functional at two different levels. The first one is the static XC potential needed 
to calculate the ground-state density of the system. The second approximation is needed to 
represent the XC kernel fxc{^, i"'; ^) which determines the XC contribution to the screening 



of an applied field. For 
approximation (ALDA) 



jhe XC kernel, we use reasonably accurate adiabatic local density 
. On the other hand, for static XC energy functional we use the 



Dirac exchange energy functional [39| and Gunnarsson-Lundqvist (GL) [40| parametrized 
form for the correlation energy functional within the LDA. In addition to approximating 
the XC energy functional, the calculation of polarizability by the hydrodynamical approach 
of TDDFT also requires approximating the non-interacting kinetic energy functional Tg[p\ 
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also. We choose the von Weizsacker 4l| form for Ts[p], which is given as 



Tw[p] = g y dr. (5) 

The von Weiszacker functional is well suited for the description of response properties of 



closed shell metal clusters 

In the present paper the ground-state densities p^^\r) of clusters are obtained by em- 
ploying purely density-based ETF jlo|, Q method within the S JBM of metal clusters. 



B. Ab initio TDDFT based method 



To carry out the all-electron ab initio calculation of the frequency dependent polariz- 



abilities of the clusters, we employ the ADF program package 421]. We refer the reader to 



Ref. [37| for detailed description of the method adopted in this package for obtaining the 
frequency dependent polarizabilities. In this package calculations of electronic and response 
properties of molecules are carried by using Slater type orbital (STO) basis sets. It is well 
known that for accurate calculations of response properties it is necessary to have large basis 
sets with both polarization and diffuse functions. For our purpose, we choose all electron 
even tempered basis set ET-QZ3P-2DIFFUSE with two sets of diffuse functions consisting 
of (lls,9p,7d,3f) functions for Na atom. The application of basis set with diffuse functions 
often leads to the problem of linear dependencies. Such problem have been circumvented by 
removing linear combinations of functions corresponding to small eigenvalues of the overlap 
matrix. We expect that the size of the chosen basis set will make our results very close to 
the basis-set limit. Once the average polarizability at imaginary frequencies are obtained 
then the Casimir-Polder integral Eq. ([2]) is evaluated by employing thirty point Gauss- 



Chebyshev quadrature scheme as described in Ref. [43|. The convergence of the results has 
been checked by increasing number of frequency points. Like SJBM based calculations, all 
the ab initio calculations in the present paper are performed with the ALDA for the XC 
kernel. On the other hand, for the static XC potential needed to calculate the ground-state 
orbitals and their energies, we employ the LDA potential as parametrized by Vosko, Wilk 



and Nussair (VWN) [4J]. We note here that both GL and VWN forms for the XC potential 



use the same Dirac exchange energy functional 39(] but the parameterization for the correla- 



tion part is different. We expect that this difference in the correlation energy functional will 
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lead to a significantly smaller deviation in the results obtained via the ab initio and SJBM 
calculations than the discrepancy arising due to the consideration of the actual geometrical 
structure of the clusters in the ab initio calculations. Furthermore, in order to perform ab 
initio calculation of response properties we need to choose the ground-state geometries of 
clusters. For the dimer Na2, we use the experimental bond length 3.0786 A. On the other 
hand, for larger clusters (4- to 20-atom clusters) we use the structures which are obtained 
via geometry optimization calculation with triple-fplus two polarization functions (TZ2P) 
basis set and Becke-Perdew (BP86) XC potential 45|, |46|] . All the optimizations are carried 
out with the convergence criteria for the norm of energy gradient and energy fixed at 10~^a.u 
and 10~^a.u, respectively. The optimized structures obtained by us are in agreement with 



the corresponding results of Refs. 



25 



261] . In case of a cluster having more than one isomers 



we choose the one possessing the lowest energy for our calculations of the dipole polarizabil- 
ity. Here we note that our geometry optimization calculation for the cluster containing 20 
sodium atoms yields lower energy for the structure with C2v symmetry than the one with 
Td symmetry. This is in contrast to the result of Ref . 25|] . The next section is devoted to 
the discussion of results for the dispersion coefficient Cq obtained by employing the above 
two methods. 



III. RESULTS AND DISCUSSION 

We begin this section with a comparison of the results for Cg between similar pairs of 
sodium clusters (Na„-Na„) with even number of atoms n ranging from 2 to 20. These results 
are shown in Table I. For comparison, we also include in Table I the results for Cr of 2-, 
8-, and 20-atom clusters obtained by employing TDKS formalism within the SJBM |7|. We 
note that the values of Cq for the dimer obtained within SJBM by employing hydrodynamic 
and TDKS approaches are almost the same. This is because the von Weizsacker form for 
the kinetic energy functional is exact for two-electron systems. For 8- and 20-atom clusters 
the difference^ et ween the TDKS and hydrodynamic results are of the order of 10%. In our 
earlier study [9] , we have found similar order of difference between the two results for 40-atom 
cluster also. This then demonstrates that the SJBM-hydrodynamic results are reasonably 
accurate and capable of yielding Cq quite close to the more accurate orbital based TDKS 
number. Having assessed the accuracy of the SJBM-hydrodynamic results for the van der 
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TABLE I: Comparison of dispersion coefficient Cg (in atomic units) for pairs of similar sodium 
clusters obtained with the all-electron ab initio and SJBM calculations. 



N 


Ce X 10 ^ a.u. 




ab initio 


SJBM-Hydrodynamic 


SJBM-KS" 


2 


3.68 


2.63 


2.62 


4 


15.05 


9.73 


- 


6 


29.44 


20.65 


- 


8 


41.82 


35.43 


40.06 


10 


71.11 


54.04 




12 


103.74 


76.43 




14 


137.45 


102.56 




16 


160.38 


132.38 




18 


191.23 


165.87 




20 


232.81 


202.99 


228.58 



(a) Ref. 



Waals coefficient Ce with respect to the corresponding TDKS results, we next focus our 
attention on their comparison with the corresponding ab initio numbers. First it is evident 
from Table I that the results based on SJBM are lower than the ab initio numbers for all the 
clusters considered in this paper. The deviation is more for the non-magic number clusters. 
In order to estimate the deviation between the two results, we plot in Fig. 1 the relative error 
given by (Cg** — Cq'^^^^)/Cq'' as a function of the number of atoms. Fig. 1 clearly 
reveals some important results of this paper. We notice that the relative error shows an 
oscillatory behaviour with the increasing number of atoms. It is important to note here that 
the relative errors for the closed shell magic number clusters assume the positions of minima. 
Moreover, the relative error for these magic number clusters decrease as the number of atoms 
present in the cluster increases. For example, the relative error for the dimer is about 28% 
and it goes down to approximately 12% for the 20-atom cluster. This is an encouraging 
results as for the clusters larger than Na2o performing ab initio calculations of response 
properties hke frequency dependent polarizability is computationally quite demanding. For 
other clusters the relative error in magnitude with maxima at 4- and 12-atom clusters. 
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"lO I I I I I I I I I I I I I I I I I I I I I I 

2 4 6 8 10 12 14 16 18 20 22 

Number of Atoms 



FIG. 1: Percent difference 100 x {Cf - Cl'^^^)/Cf between the ah initio and SJBM- 
hydrodynamic results for the van der Waals coefficients as a function of number of atoms. The 
lines joining the points are guide to the eye. 

However, the value of relative error for the 12-atom cluster is lower than for the 4-atom one. 
We expect that for larger non-magic number clusters the SJBM results for Ce will be closer 
to the corresponding ah initio numbers. Encouraged by these results, we perform calculation 
of Ce for all pairs of clusters {Nan — Nan). The results of these calculations are presented 
in Table II. Again it is seen from Tabic II that for different pairs of magic-number clusters 
the relative difference between the ah initio and the SJBM results is less as compared to the 
other pairs. From Table II, we also infer that the jellium based results for clusters larger 
than the 14-atom one are quite close to the corresponding ah initio numbers. These results 
then demonstrate that for clusters larger than Nai4 the detailed ionic core structure does 
not have much influence on the values of the van der Waals coefficient Ce. 
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TABLE II: Comparison of dispersion coefficient Cg xlO~^ (in atomic units) for all pairs of sodium 
clusters obtained with the all-electron ab initio and SJBM calculations. The SJBM results are 
shown in the parenthesis 



N 


2 


4 


6 


8 


10 


12 


14 


16 


18 


20 


2 


3.68 


7.45 


10.41 


12.39 


16.18 


19.54 


22.49 


24.27 


24.48 


29.22 




(2.63) 


(5.06) 


(7.38) 


(9.66) 


(11.93) 


(14.19) 


(16.43) 


(18.67 ) 


(20.89) 


(23.11) 


4 




15.05 


21,04 


25.02 


32.67 


39.47 


45.42 


48.99 


53.44 


63.85 






(9.73) 


(14.17) 


(18.56) 


(22.92) 


(27.25) 


(31.57) 


(35.87 ) 


(40.15) 


(44.41) 


6 






29.44 


35.05 


45.74 


55.25 


63.59 


68.59 


74.91 


89.40 








(20.65) 


(27.05) 


(33.41) 


(39.73) 


(46.02) 


(52.28) 


(58.52) 


(64.74) 


8 








41.82 


54.51 


65.84 


75.80 


81.89 


89.42 


98.66 










(35.43) 


(43.76) 


(52.04) 


(60.28) 


(68.48) 


(76.66) 


(84.75) 


10 










71.11 


85.89 


98.86 


106.75 


116.53 


128.57 












(54.04) 


(66.26) 


(74.44) 


(84.57) 


(94.67) 


(104.73) 


12 












103.74 
(76.43) 


119.41 
(88.53) 


128.93 
(100.59) 


140.74 
(112.59) 


155.28 
(124.55) 


14 














137.45 
(102.56) 


148.43 

(116.52) 


162.03 

(130.43) 


178.77 

(144.28) 


16 
















160.38 
(132.38)) 


175.12 
(148.18) 


193.22 
(163.93) 


18 


















191.23 
(165.87) 


211.06 

(183.50) 


20 




















232.81 
(202.99) 



Having assessed the accuracy of results obtained within the SJBM, we now focus our 
attention on the study of relationship between Cg and Oi{0). In the so-called London approx- 
imation the van der Waals coefficient Cg between two similar molecules can be represented 
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in terms of the static polarizability a(0) as 

a = ^(a(0)f, (6) 

where ui is the effective or characteristic frequency of the system. The above expression (Eq. 
([6])) is obtained with the single pole approximation for the frequency dependent polarizabil- 
ity, which assumes that all the transitions can be replaced by only one effective transition. 
This is also known as London dispersion formula and it provides a way to correlate the van 
der Waals coefficient with the static polarizability. In doing so the accuracy of Ce crucially 
depends on the precision with which the static polarizability is calculated. In order to study 
this correlation, we calculate the average static polarizability of sodium clusters by the ab 
initio method and SJBM-hydrodynamic approach. Another reason for carrying out static 
polarizability calculations is that unlike Ce, a fairly large number of experimental results on 



2400 




FIG. 2: Plot of average static polarizability a(0) of sodium atom clusters in atomic units. The solid 
squares represent the ab initio results, SJBM-hydrodynamic results are denoted by solid circles, 
and solid triangles represent the experimental results of Ref. [47i]. The lines joining the points are 
guide to the eye. 
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the static polarizabilty exist in the hterature 47|, |48|], which allows us to compare our results 
with the corresponding experimental data. The results of the polarizability calculations are 
shown in Fig. 2 along with the experimental results of Ref. 47|. Fig. 2 clearly shows 
that both ab initio and SJBM results are lower than the corresponding experimental data, 
however, the ab initio results are closer to the experimental data than the SJBM results. 
Moreover, it is seen from this figure that like Cq, the deviation between the SJBM and ab 
initio results is minimum for magic number clusters and these deviations decrease with the 
larger clusters. For example, the difference between the two results is 20% ( ab initio number 
is 233.32 a.u.) for the dimer, whereas the differnce reduces to around 11% ( ab initio number 
is 1752.8 a.u.) for the Na2o cluster. For non-magic number clusters the differnce is more and 
but reduces with the increase in the size of the clusters. Therefore, we conclude that the 
static polarizability of large size clusters can also be obtained quite accurately by employing 
the SJBM-hydrodynamic approach with substantially less amount of computational effort. 



IV. CONCLUSION 



The van der Waals coefficient Ce for sodium atom clusters containing even number of 
atoms ranging from 2 to 20 atoms have been calculated by employing all-electron ab initio 
method and the SJBM-hydrodynamic formalism within the realm of TDDFT. The calcula- 
tions are performed with the LDA XC potentials. The van der Waals coefficient is obtained 
by using Casimir-Polder expression which requires frequency dependent dipole polarizabil- 
ties of the two interacting species at imaginary frequencies. These are calculated by the 
all-electron ab initio and the SJBM-hydrodynamic methods of TDDFT. The results of these 
two calculations are compared to ascertain the effect of ionic structure of the clusters on the 
van der Waals coefficient. We find that the differences between the results of all-electron ab 
initio and SJBM-hydrodynamic calculations are minimum for magic number clusters and 
these differences decrease with the increase in the number of atoms present in the clusters. 
Even for non-magic number clusters, we find that SJBM-hydrodynamic calculations lead to 
results which are very close to the corresponding ab initio numbers. From these results, we 
conclude that for clusters larger than Nai4 detailed ionic structure of the clusters does not 
have much effect on the results for the van der Waals coefficient and also on their static 
polarizability. This is an important results as carrying out ab initio calculations for the 
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optical response properties of large clusters is very expensive. This paper shows that for 
such systems SJBM based methods within TDDFT can yield quite accurate results for the 
van der Waals coefficient. 
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